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ABSTRACT 


The  development  of  accurate  numerical  methods  to  simulate  wave  propagation  in  three-dimensional  (3D)  earth 
models  and  advances  in  computational  power  offer  exciting  possibilities  for  modeling  the  motions  excited  by 
underground  nuclear  explosions.  This  presentation  will  describe  recent  work  to  use  new  numerical  techniques  and 
parallel  computing  to  model  earthquakes  and  underground  explosions  to  improve  understanding  of  the  wave 
excitation  at  the  source  and  path-propagation  effects.  Firstly,  we  are  using  the  spectral  element  method  (SEM, 
SPECFEM3D  code  of  Komatitsch  and  Tromp,  2002)  to  model  earthquakes  and  explosions  at  regional  distances 
using  available  3D  models.  SPECFEM3D  simulates  anelastic  wave  propagation  in  fully  3D  earth  models  in 
spherical  geometry  with  the  ability  to  account  for  free  surface  topography,  anisotropy,  ellipticity,  rotation  and 
gravity.  Results  show  in  many  cases  that  3D  models  are  able  to  reproduce  features  of  the  observed  seismograms  that 
arise  from  path-propagation  effects  (e.g.,  enhanced  surface  wave  dispersion,  refraction,  amplitude  variations  from 
focusing  and  defocusing,  tangential  component  energy  from  isotropic  sources).  We  are  currently  investigating  the 
ability  of  different  3D  models  to  predict  path-specific  seismograms  as  a  function  of  frequency.  A  number  of  models 
developed  using  a  variety  of  methodologies  are  available  for  testing.  These  include  the  WENA/Unified  model  of 
Eurasia  (e.g.  Pasyanos  et  al  2004),  the  global  CUB  2.0  model  (Shapiro  and  Ritzwoller,  2002),  the  partitioned 
waveform  model  for  the  Mediterranean  (van  der  Lee  et  al.,  2007)  and  stochastic  models  of  the  Yellow  Sea  Korean 
Peninsula  region  (Pasyanos  et  al.,  2006).  Secondly,  we  are  extending  our  Cartesian  anelastic  finite  difference  code 
(WPP  of  Nilsson  et  al.,  2007)  to  model  the  effects  of  free-surface  topography.  WPP  models  anelastic  wave 
propagation  in  fully  3D  earth  models  using  mesh  refinement  to  increase  computational  speed  and  improve  memory 
efficiency.  Thirdly,  we  are  modeling  non-linear  near-source  shock  wave  propagation  with  GEODYN,  an  Eulerian 
Godunov  finite-difference  code  (Antoun  et  al.,  2001).  This  code  accounts  for  shock  wave  propagation  and  a  variety 
of  effects,  including  cavity  formation,  rock  fracture  and  plastic  deformation.  We  are  exploring  the  coupling  of 
GEODYN  to  WPP  to  propagate  motions  from  the  near-source  (non-linear)  region  to  the  (linear  anelastic)  region 
where  seismic  observations  are  made  at  local,  regional  and  teleseismic  distances.  This  effort  has  just  begun  and  we 
show  preliminary  results  in  this  paper.  These  simulation  tools  are  supported  by  massively  parallel  computers 
operated  by  Livermore  Computing. 
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OBJECTIVES 

Improvements  in  numerical  methods  to  simulate  wave  phenomena  in  three-dimensional  (3D)  heterogeneous 
geologic  materials,  implementation  of  these  algorithms  on  parallel  computers  and  the  inexorable  increase  of 
computational  power  offer  exciting  possibilities  for  the  ability  to  simulate  wave  propagation  for  improving 
understanding  of  the  nuclear  explosion  source  and  nuclear  explosion  monitoring.  The  objective  of  this  study  is  to 
describe  the  development  of  new  seismic  wave  modeling  capabilities  at  Lawrence  Livermore  National  Laboratory 
(LLNL)  and  articulate  how  these  capabilities  can  improve  understanding  of  observed  seismic  waves,  particularly 
source  excitation  and  propagation  effects.  Nuclear  explosion  monitoring  based  on  seismic  data  ultimately  relies  on 
understanding  the  fundamentals  of  excitation  of  waves  in  the  earth  by  the  high-energy  nuclear  explosion  source, 
shock-wave  propagation  in  the  near-source  region  and  subsequent  elastic  wave  propagation,  including  interaction 
with  near-source  geology  and  topography  and  path  propagation.  We  are  applying  newly  developed  numerical 
methods  and  high  performance  computational  resources  to  improve  knowledge  of  observable  effects  on 
seismograms  due  to  propagation  in  realistic  3D  earth  models  and  geologic  materials. 


RESEARCH  ACCOMPLISHED 

This  effort  covers  three  areas  of  research:  1)  simulation  of  regional  distance  (<  1500  km)  seismograms  at 
intermediate  periods  (0.02-0.1  Hz)  with  the  spectral  element  method  (SEM);  2)  simulation  of  realistic  free-surface 
topography  with  our  Cartesian  anelastic  finite  difference  code  (WPP);  and  3)  the  simulation  of  non-linear  shock 
wave  phenomena  in  realistic  geologic  materials  (GEODYN)  and  the  coupling  of  non-linear  motions  to  our  linear 
anelastic  code  (WPP).  All  these  efforts  rely  on  massively  parallel  computer  systems  at  LLNL. 

Simulations  of  Seismic  Wave  Propagation  in  3D  Models  at  Regional  Distances  and  Intermediate  Periods 

Researchers  in  the  academic  and  nuclear  monitoring  communities  have  vigorously  developed  3D  seismic  velocity 
models  based  on  various  seismic  tomographic  methods.  Computational  methods  for  seismic  wave  propagation  using 
parallel  computers  allow  us  to  compute  synthetic  seismograms  in  3D  models  for  intermediate  periods  (>  10  seconds) 
at  regional  distances  (<  1500  km).  To  do  this,  we  are  using  the  spectral  element  method  (SEM)  code  SPECFEM3D 
(Komatitsch  and  Vilotte,  1998;  Komatitsch  and  Tromp,  1999;  Komatitsch  et  al.,  2000).  This  code  computes  the 
response  of  the  earth  to  arbitrary  forcing  in  spherical  geometry  with  fully  3D  seismic  velocity  (including  anisotropy) 
and  density  variations.  The  code  also  handles  the  effects  of  free-surface  topography  and  bathymetry,  ellipticity, 
rotation  and  self- gravitation.  We  have  modified  the  code  to  handle  a  simple,  but  general  3D  model  format. 


Figure  1.  (left)  SPECFEM3D  computational  domain  for  regional  simulations  of  events  in  East  Asia,  (right) 
Map  showing  the  locations  of  the  January  11,  2000  (2000/011)  earthquake  and  October  9,  2006, 
(2006/282)  North  Korean  nuclear  explosion  along  with  regional  broadband  stations  that  recorded 
the  events. 
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We  have  used  this  code  to  compute  synthetics  seismograms  for  the  January  11,  2000,  moderate  (Mw  ~  5)  earthquake 
in  the  Korean  Peninsula  and  the  (Mw  -  4)  October  9,  2006,  North  Korean  nuclear  explosion.  Figure  1  (left)  shows 
SPECFEM3D  computational  domain  at  the  surface  (white  box)  covering  a  45°  “chunk”  (solid  angle)  of  the  earth. 
Figure  1  (right)  shows  a  map  of  the  events  and  recording  stations  for  these  two  events.  Focal  parameters  for  the 
earthquake  were  obtained  from  waveform  modeling  using  ID  methods. 


Time,  seconds  Time,  seconds 

Figure  2.  Observed  (blue)  and  synthetic  seismograms  for  the  2000/011  earthquake  at  station  BJT  for  the 
iasp91  (green)  and  s29ea+CRUST2.0  (red)  models  filtered  0.02-0.05  Hz  (left)  and  0.03-0.1  Hz 
(right). 

Figure  2  shows  the  observed  and  synthetic  waveforms  for  the  2000/01 1  earthquake  at  station  BJT  (Beijing,  China) 
filtered  in  two  different  bands:  0.02-0.05  and  0.03-0.07  Hz  (left  and  right,  respectively).  Synthetics  were  computed 
for  two  models:  the  iasp91  ID  model  (Kennett  and  Engdahl,  1991)  shown  as  the  green  synthetics  and  the 
s29ea+CRUST2.0  3D  model  (Kustowski  et  al.,  2008)  shown  as  the  red  synthetics.  Calculations  resulted  in 
seismograms  that  were  resolved  in  the  frequency  band  0.0-0.12  Hz  and  took  about  4  hours  on  64  CPU’s  of  a  modest 
Linux  cluster  operated  by  the  Ground-Based  Nuclear  Explosion  Monitoring  (GNEM)  Program  at  LLNL.  All 
calculations  included  the  effects  of  surface  topography,  bathymetry  and  ellipticity.  This  event  was  moderately  large 
(Mw  ~  5.0)  and  had  good  signal-noise  across  a  broadband.  Note  that  for  0.02-0.05  Hz  (20-50  second  periods)  the 
observed  and  synthetic  waveforms  are  very  simple  with  clear  Pnl,  Rayleigh  and  Love  waves.  At  these  periods  the 
synthetics  match  the  observed  waveforms  with  the  iasp9 1  model  showing  a  slightly  better  fit  to  the  phase  and 
amplitudes  of  the  data.  However,  including  slightly  higher  frequencies  (0.1  Hz  or  10  seconds  periods)  the  observed 
waveforms  show  much  greater  complexity,  especially  the  surface  waves.  This  complexity  is  likely  due  to  interaction 
with  the  Bohai  Basin  (Figure  1).  Note  that  the  Rayleigh  wave  is  strongly  dispersed  and  the  later  part  of  the  Love 
wave  shows  a  long  duration  with  little  dispersion.  The  s29ea+CRUST2.0  3D  model  captures  some  of  the  more 
complex  surface  waveform,  mostly  in  the  Rayleigh  wave.  As  expected  at  intermediate  frequencies  (0.02-0.07  Hz  or 
14-20  second  periods)  the  observed  waveforms  are  less  complex. 

This  example  (Figure  2)  demonstrates  the  challenge  for  waveform-based  monitoring  of  small-to-moderate  sized 
events  (Mw  3. 5-4. 5).  These  events  radiate  very  little  long-period  energy  and  are  often  observed  where  the  surface 
waves  only  have  signals  above  the  noise  for  intermediate  periods  (e.g.,  0.03-0.1  Hz  or  10-33  seconds).  For 
waveform  analysis  at  intermediate  periods  seismologists  require  good  constraints  on  the  path-specific  structure.  We 
are  performing  calculations  with  a  variety  of  available  3D  models  to  determine  the  performance  of  models  for 
predicting  observed  waveforms  as  a  function  of  frequency  content.  It  is  important  to  quantify  the  performance  of  3D 
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models  so  that  they  may  be  used  as  starting  models  for  waveform-based  tomography,  such  as  adjoint  methods 
(Tromp  et  ah,  2005;  Liu  and  Tromp,  2006). 
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Figure  3.  Observed  (black)  and  synthetic  (red)  three-component  seismograms  for  the  ID  iasp91  (left)  and  3D 
CUB2.0  (right)  models.  Both  data  and  synthetics  are  filtered  0.05-0.1  Hz  (10-20  second  periods). 


As  another  example  of  the  effects  of  known  3D  structure  on  intermediate  period  surface  waves  we  show  the 
observed  and  synthetic  waveforms  for  the  October  9,  2006,  North  Korean  nuclear  test  at  station  BJT.  The  path  to 
BJT  passes  through  the  northern  edge  of  the  Bohai  Basin,  similar  to  the  path  for  the  2000/01 1  earthquake  (Figure  1). 
The  synthetics  for  the  ID  iasp91  and  3D  CUB2.0  model  (Shapiro  and  Ritzwoller,  2002)  are  shown  in  Figure  3. 
While  the  signal-to-noise  ratio  of  the  transverse  component  at  the  expected  time  window  of  the  surface  wave  is 
about  1:1  or  slightly  higher,  the  3D  synthetic  predicts  some  energy  that  is  not  present  in  the  synthetics  for  a  ID 
model.  This  illustrates  how  the  refraction  of  surface  waves  by  3D  structure,  such  as  sedimentary  basins,  can  cause 
rotation  of  energy  from  a  purely  isotropic  source  from  the  radial  to  the  transverse  component.  Such  wave 
propagation  effects  may  bias  estimates  of  source  parameters  for  explosion  sources  such  as  seismic  moment  (and 
consequently  explosive  yield)  and  deviatoric  moment  tensor  components. 
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Figure  4.  WENA  3D  model  properties  across  the  Middle  East:  sediment  thickness  (top  left)  and  crustal 
thickness  (bottom  left).  Snapshots  of  vertical  ground  displacements  showing  Rayleigh  wave 
propagation  for  an  isotropic  source  using  an  average  ID  model  (center)  and  the  WENA  3D  model 
(right). 
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Finally,  we  performed  large-scale  simulations  of  regional  seismic  wave  propagation  in  the  Middle  East  using  our 
anelastic  finite  difference  code  (WPP,  described  below).  We  modeled  the  Mw  6.25,  February  22,  2005,  Kerman 
earthquake  using  an  average  ID  model  and  the  3D  WENA  model  of  Pasyanos  et  al.  (2004).  The  WENA  model 
represents  3D  crustal  and  upper  mantle  structure  with  1°  resolution  and  reveals  the  known  sedimentary  and  crustal 
structure  of  the  Middle  East  (Figure  4,  left).  For  simplicity  of  illustration,  we  compare  the  wave  propagation  excited 
by  an  isotropic  source  using  the  ID  and  WENA  3D  models  (Figure  4,  center  and  right,  respectively).  The  vertical 
displacement  wavefields  are  shown  at  two  different  time  steps  in  the  calculations.  Notice  the  dramatic  timing, 
amplitude  and  dispersion  asymmetries  in  the  3D  model  compared  to  the  ID  model.  These  effects  are  confirmed  by 
observations  of  long-duration  intermediate  period  surface  waveforms  for  events  in  the  Zagros  Mountains  and  Iranian 
Plateau  recorded  in  the  eastern  Arabian  Peninsula  along  the  Persian/ Arabian  Gulf  coast. 

Note  that  the  simulations  shown  in  Figure  4  required  nearly  4  billion  grid  points  for  density,  seismic  velocities  and 
quality  factors  covering  a  volume  3,000  km  by  3,000  km  by  400  km  (depth)  at  1,000-m  resolution.  The  calculations 
took  nearly  50,000  CPU-hours,  taking  nearly  two  days  on  1,024  CPU’s  of  the  Thunder  Linux  cluster  at  LLNL.  The 
simulations  resulted  in  seismograms  valid  to  0.1  Hz  (10  second  periods). 

Simulation  of  Seismic  Wave  Propagation  with  Realistic  Free  Surface  Topography 

LLNL  recently  developed  a  new  parallel  finite  difference  code  for  simulating  anelastic  wave  propagation.  The 
algorithm  uses  a  node-center  second-order  discretization,  described  by  Nilsson  et  al.  (2007).  The  code,  called  WPP, 
is  open  source  and  available  at  https://computation.llnl.gov/casc/serpentine/index.html.  Currently,  the  code  handles 
fully  3D  variations  in  seismic  wave  speeds  and  density,  both  acoustic  and  elastic  wave  propagation  and  attenuation 
with  quality  factors  for  P-  and  S-waves.  The  code  allows  for  mesh  refinement  for  memory  and  computational 
efficiency.  As  the  resolvable  frequency  of  seismic  events  increases,  the  corresponding  wavelength  of  the  seismic 
wave  motion  decreases  requiring  finer  discretization  and  more  grid  points  (and  memory).  For  frequencies  above 
about  0.5  Hz,  effects  of  non-planar  topography  become  increasingly  important  (Figure  5).  We  are  currently 
extending  the  WPP  code  to  accurately  satisfy  the  free  surface  boundary  conditions  on  realistic  non-planar 
topography.  This  will  be  accomplished  by  using  an  energy  conserving,  summation-by-parts,  finite  difference 
discretization  of  the  elastic  wave  equation  on  a  curvilinear  grid  (Appelo  and  Petersson,  2008). 


Figure  5.  (left)  Vertical  displacement  time  history  along  the  surface  due  to  an  isotropic  point  source  below  the 
surface.  Note  the  effects  of  topography  on  the  right  side  of  the  domain,  (right)  Compressional 
velocity  (VP)  in  an  east-west  cross-section  of  the  San  Francisco  Bay  area  going  through  Mt. 
Tamalpais.  The  scale  goes  from  5,000  m/s  (dark  red)  to  0  m/s  (dark  blue).  The  cross  section  extends 
to  a  depth  of  10  km  below  sea  level. 

Since  the  curvilinear  formulation  involves  about  nine  times  more  arithmetic  operations  per  grid  point  than  its 
Cartesian  counterpart,  the  computational  grid  will  be  made  Cartesian  below  some  fixed  elevation.  As  is  visible  in 
Figure  5  (right),  the  compressional  wave  speed  (VP)  increases  drastically  with  depth  and  is  commonly  about  eight 
times  larger  below  the  Moho  discontinuity  (at  about  30  km  depth)  than  in  sedimentary  basins  near  the  free  surface. 
The  propagation  speed  for  shear  waves  has  similar  depth  dependence,  but  the  ratio  between  the  two  speeds  depends 
on  the  material  type  (sediment,  hard  rock,  etc.).  Since  the  wavelength  of  elastic  waves  is  proportional  to  the  wave 
speed,  structured  mesh  refinement  (where  the  mesh  gets  coarser  at  depth)  will  be  used  to  keep  the  grid  size  in 
approximate  parity  with  the  wavelength.  Furthermore,  since  there  is  such  a  large  change  in  wave  speeds  from  the 
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surface  to  the  bottom  of  the  computational  domain,  mesh  refinement  will  significantly  reduce  the  total  number  of 
grid  points  and  also  allow  the  explicit  time-stepping  to  use  a  much  larger  time  step  (improving  computational 
efficiency).  To  guarantee  stability  of  the  time-stepping  scheme  in  the  presence  of  mesh  refinement,  we  will  use  a 
newly  developed  energy  conserving  coupling  of  the  solution  across  mesh  refinement  boundaries.  For  simplicity, 
only  factor-of-two  refinements  are  allowed  and  all  mesh  refinement  boundaries  must  be  horizontal.  The  latter 
restriction  allows  for  a  two-dimensional  decomposition  of  the  computational  domain  where  each  processor  has  a 
“pencif’-shaped  domain  from  the  topography  at  the  top  to  the  far- field  boundary  at  the  bottom.  This  layout  allows 
for  efficient  load  balancing  on  parallel  machines  and  eliminates  the  need  for  expensive  all-to-all  communications. 

The  influence  of  topography  is  illustrated  in  Figure  6,  where  we  show  the  response  due  to  a  vertically  propagating 
shear  wave  impinging  on  a  3-D  Gaussian  hill.  In  this  case,  the  material  is  homogeneous,  a  free  surface  boundary 
condition  is  imposed  on  the  curved  top  surface,  and  periodic  boundary  conditions  are  used  in  the  horizontal 
directions.  The  traces  show  displacement  time  histories  at  the  surface  along  a  densely  sampled  line  passing  through 
the  center  of  the  Gaussian  hill.  For  a  flat  free  surface,  the  horizontal  component  would  only  have  one  pulse 
corresponding  to  the  arrival  of  the  vertically  propagation  wave;  the  vertical  component  would  be  identically  zero.  In 
comparison,  the  curved  topography  induces  reflected  P-  and  S-waves,  which  also  induce  a  vertical  component  to  the 
motion. 
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Figure  6.  Horizontal  (left)  and  vertical  (right)  displacement  time  histories  along  a  line  through  the  center  of 
the  Gaussian  hill. 


Simulation  of  Non-Linear  Hydrodynamics  with  GEODYN  and  Coupling  to  WPP 

The  high  energy  density  of  chemical  and  nuclear  explosions  causes  irreversible  processes  (e.g.,  vaporization,  cavity 
creation  and  plastic  deformation)  in  geologic  materials  near  the  source  before  stimulating  seismic  waves  beyond  the 
elastic  radius.  In  order  to  simulate  ground  motions  emerging  from  buried  explosions,  we  must  account  for  these 
processes  including  a  complete  thermodynmical  treatment  of  the  energy  involved  in  these  processes.  Source 
emplacement  conditions  are  known  to  have  dramatic  effects  on  the  excitation  of  seismic  waves  that  ultimately  must 
be  accounted  for  in  yield  estimates  from  seismic  amplitudes  and  earthquake-explosion  discriminants.  The  challenges 
of  site-specific  yield  scaling  relationships  and  the  transportability  of  high-frequency  regional  discriminants  can  be 
traced  to  emplacement  conditions  and  physical  properties  of  geologic  materials  at  the  source. 

The  subject  of  this  effort  is  to  develop  and  illustrate  a  coupling  approach  between  the  non-linear  GEODYN 
hydrodynamics  code  (Antoun  et  al.,  2001)  and  the  linear  anelastic  wave  propagation  code  WPP,  described  above 
(Nilsson  et  al.,  2007).  GEODYN  is  a  nonlinear  Eulerian  wave  propagation  code  with  adaptive  mesh  refinement 
(Antoun  et  al.,  2001;  Lomov  and  Rubin,  2003).  The  GEODYN  code  can  handle  large  material  deformations  and 
phase  transitions  such  as  melting,  vaporization,  porous  compaction  and  material  damage.  In  the  limit  of  low 
compressions  and  temperatures  GEODYN  can  be  applied  to  model  linear  waves  in  solid  media,  albeit  at  significant 
numerical  cost  over  conventional  anelastic  finite  difference  methods.  The  code  is  not  very  efficient  for  handling 
linear  waves  since  it  carries  around  a  lot  of  extra  variables  and  performs  calculations  that  are  not  needed  for  linear 
(elastic)  wave  propagation.  WPP,  on  the  other  hand,  cannot  handle  permanent  deformations  and  uses  a  simple 
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scheme  to  model  inelastic  effects  on  the  large  scale  (Liu  and  Archuleta,  2006).  The  purpose  of  this  study  (which  has 
only  just  begun)  is  to  couple  the  two  codes  to  model  the  non-linear  behavior  of  the  near-source  region  with  the  full 
physics  capabilities  of  GEODYN  and  pass  these  motions  to  WPP  in  an  overlapping  region  in  order  to  propagate  the 
seismic  waves  to  local,  regional,  and/or  teleseismic  distances.  The  basic  idea  of  coupling  is  shown  schematically  in 
Figure  7. 


Figure  7.  Conceptual  illustration  of  the  coupling  between  the  nonlinear  source  region  from  GEODYN  (GD)  to 
WPP.  The  GEODYN  domain  extends  to  radius  R  (dark  gray  region),  while  the  WPP  domain  covers 
a  larger  entire  volume. 

GEODYN  can  compute  motions  in  ID  (spherically  symmetric),  2D,  2.5D  (2D  axisymmetric)  or  fully  3D  domains. 
Figure  7  describes  a  spherical  source  modeled  in  a  ID  GEODYN  calculation  that  was  passed  to  a  3D  WPP 
calculation.  Material  displacement  inside  the  volume  within  the  hand-off  radius  R  was  overwritten  with 
pre-calculated  ID  GEODYN  data.  One  of  the  most  important  applications  for  GEODYN  is  the  accurate  simulation 
of  ground  motion  due  to  underground  explosions.  Experimental  data  on  spherical  waves  generated  in  situ  by 
underground  explosions  published  in  open  literature  are  scarce  and  are  not  always  available  with  good  geotechnical 
characterization  of  the  emplacement  rock.  For  this  study,  we  chose  to  use  data  from  a  series  of  chemical 
high-explosive  (HE)  experiments  in  limestone  performed  in  Kirghizia  in  1960  (Murphy  et  al.,  1997)  to  illustrate  the 
code-coupling  approach.  The  generic  strength  model  described  in  Vorobiev  (2008)  was  used  in  GEODYN  to  model 
0.5%  porous  limestone.  Experimental  measurements  of  ground  motion  at  different  distances  from  the  shot  point 
were  reported  by  Murphy  et  al.  (1997).  The  explosions  and  motion  recordings  reported  in  this  study  were  conducted 
in  the  subsurface  in  essentially  rock  whole-space  conditions,  making  them  ideal  for  our  simulation  experiments. 

Figure  8  shows  the  vertical  component  displacement  time  histories  calculated  in  ID  GEODYN  (blue)  and  passed  to 
WPP  (red)  at  the  boundary  of  the  hand-off  region  (point  B  in  Figure  7).  The  spatial  resolution  used  in  GEODYN 
was  6  cm  and  the  space  resolution  used  in  3D  WPP  was  1  m,  allowing  for  resolution  to  about  100  Hz.  Elastic 
constants  used  in  WPP  were  chosen  to  match  initial  moduli  for  the  generic  limestone.  Note  that  the  amplitude  and 
phase  agree  very  well. 

Figure  9  shows  peak  displacements  as  a  function  of  range  calculated  by  ID  GEODYN  and  3D  WPP  codes.  The  ID 
GEODYN  calculations  agree  very  well  with  experimentally  measured  values  of  peak  velocity  (crosses  for  point 
measurements  taken  from  Figure  8  of  Murphy  et  al.,  1997).  The  peak  displacements  computed  with  WPP  agree  well 
with  the  GEODYN  calculations  when  a  fine  grid  spacing  of  1  meter  is  used.  However,  the  WPP  peak  displacements 
attenuate  faster  than  1/R  when  a  grid  spacing  of  4  meters  is  used.  This  demonstrates  that  one  must  rigorously 
evaluate  the  numerical  resolution  of  such  calculations  to  be  confident  that  numerical  dispersion  is  minimized. 
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Figure  8.  Vertical  component  displacement  time  history  at  the  boundary  of  the  hand-off  region  (point  B  in 
Figure  7)  for  the  3D  WPP  (red  solid  line)  and  ID  GEODYN  (blue  dotted  line). 


CONCLUSIONS  AND  RECOMMENDATIONS 


Recent  development  of  numerical  methods  for  simulation  of  the  3D  seismic  wave  excitation  and  propagation, 
coupled  with  ever-increasing  computational  power  are  offering  exciting  new  capabilities  for  improving  knowledge 
of  seismic  wave  phenomenology  for  nuclear  explosion  monitoring.  In  this  study,  we  demonstrated  several  efforts 
underway  at  LLNL  that  promise  to  advance  capabilities  for  various  aspects  of  seismic  wave  simulation. 

The  SEM  code,  SPECFEM3D,  models  seismic  wave  propagation  in  spherical  geometry  including  important  effects 
(e.g.,  3D  heterogeneity,  anisotropy,  free  surface  topography,  ellipticity).  This  code  allows  us  to  compute 
intermediate  period  seismograms  (periods  greater  than  about  10  seconds)  at  regional  distance  (up  to  about  1500  km). 
We  are  using  this  code  to  test  current  available  3D  seismic  models.  This  code  will  be  used  in  another  study  to 
exploit  finite  frequency  sensitivity  kernels  for  seismic  waveform  tomography  based  on  adjoint  methods  (Savage  et 
al.,  these  Proceedings). 

We  are  improving  our  Cartesian  anelastic  finite  difference  code  (WPP)  to  add  the  effects  of  free-surface  topography. 
The  method  will  rely  on  a  summation-by-parts  principle  and  has  been  proven  in  research  codes.  The  production 
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version  of  the  WPP  (version  2.0)  will  include  the  effects  of  free- surface  topography  and  other  improvements.  Note 
that  WPP  is  an  open  source  code  and  can  be  downloaded  from  LLNL. 

The  coupling  of  our  non-linear  hydrodynamics  code,  GEODYN,  with  WPP  will  allow  us  to  simulate  ground 
motions  from  underground  nuclear  explosions  including  an  accurate  representation  of  the  important  near-source 
effects  of  cavity  formation,  melting,  yielding,  plastic  deformation,  porous  compaction,  and  material  damage  all  in  a 
consistent  thermodynamic  framework.  While  many  details  remain  to  be  sorted  out,  we  have  made  good  progress  in 
the  short  time  that  we  have  been  working  on  this  effort. 


These  efforts  involve  complex  numerical  codes  running  on  massively  parallel  computers.  LLNL  has  significant 
experience  in  parallel  code  development  as  well  as  the  computational  resources  to  enable  simulations  of  relevance  to 
the  nuclear  explosion  source  phenomenology  and  monitoring  investigations. 


Figure  9.  Peak  displacement  versus  range  for  a  chemical  HE  source  in  a  limestone  whole-space: 

experimentally  measured  values  (crosses)  from  Figure  8  of  Murphy  et  al.  (1997);  calculated  with  ID 
GEODYN  (solid  green  line);  and  calculated  with  3D  WPP  (red  line  for  4-m  grid  spacing  and  black 
line  for  1-m  grid  spacing).  These  are  compared  with  and  the  expected  geometric  spreading  for  a 
spherical  wave  in  a  homogeneous  whole-space  (1/R,  fine  dashed  line). 
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